The behavior of physical systems in response to excitations can often be modeled with a linear system of equations. For example, electronic or electrical circuits or networks may be modeled by Kirchhoff's laws, which relate the currents through and voltages across the network components (such as, e.g., resistors, capacitors, transmission lines, voltage sources, etc.). Excitations of the systems may be characterized in terms of input currents and/or voltages provided by current/voltage sources and/or fixed at terminals of the network. Similarly, other physical networks, such as, e.g., pneumatic or hydraulic networks, may be modeled reasonably well by appropriate linear equations for the associated physical quantities. The linear system of equations, typically expressed as a system matrix equation, may be used to simulate the way the physical quantities associated with the system behave or evolve over time.
Physical fields, i.e., physical quantities defined at each point within a user-defined domain (e.g., a surface or volume specified geometrically, and generally corresponding to a two-, or three-dimensional region of a physical object), can likewise often be simulated using a system matrix equation. For that purpose, the fields are modeled by field-governing equations and applicable boundary conditions, which, together, describe the behavior of the field within and at the boundaries of the domain. For example, Maxwell's equations describe the behavior of an electromagnetic field, and the heat equation describes the behavior of a temperature field. The field values at the domain boundaries, as well as any source (or sink) terms in the equations, constitute the “excitations” of the system. To facilitate computational solution of the problem, the field-governing equations and boundary conditions are then discretized, resulting in a (typically large) linear system of equations. Discretization facilitates finding approximate solutions of differential or integral equations (such as, e.g., Maxwell's equations), i.e., it enables problems lacking exact mathematical (“analytical”) solutions to be approximated computationally (“numerically”). A numerical technique such as, for example, a finite difference, finite volume, and/or finite element method may be used to accomplish the discretization.
To numerically simulate the physical system (e.g., the network or field), the system matrix equation may be solved with a direct or iterative solver, depending on the size and characteristics of the linear system. (A “solver,” as the term is used herein, denotes a method for solving a system of equations, or a computer program implementing such a method, as determined by context.) Whereas a direct solver typically factorizes the system matrix into a product of special matrices (such as, e.g., the upper and lower triangular portions of the system matrix) and, thereafter, solves the equations by forward and/or backward substitution, an iterative solver solves the equations simultaneously in steps, with each step refining a previous approximation to more closely approach the exact solution. For large problems (e.g., three-dimensional electromagnetic problems defined over a large domain), direct solvers potentially require prohibitive amounts of memory and suffer poor parallel scalability. Therefore, iterative solvers typically present the only practical means for solving large systems. The iterative process is typically interrupted when a “residual vector” that quantifies the deviation of the approximation from the exact solution falls below a user-defined error, indicating the acceptability of the approximation—i.e., convergence of the solution. In order to increase the convergence rate of the system so that fewer iterations are needed to reach the acceptable approximation, the matrix may be “preconditioned” by applying a “preconditioner” matrix to the system matrix.
Many problems can be solved efficiently using a sub-class of iterative solvers usually referred to as “Krylov subspace iterative solvers.” In Krylov subspace methods, the solution is sought in a vector subspace whose dimensionality grows by one with each iteration. Conceptually, as the dimensionality of the subspace is gradually increased, the error components of the approximation correspondingly decrease. The Krylov subspace is spanned by Krylov vectors, which may be chosen to be, or based on, the residual vectors associated with each iteration. Since the residual vectors depend on the excitation, problems with multiple excitations are conventionally solved by constructing the Krylov subspace, and computing the solution vector therein, for each excitation independently. This generic approach is, however, often computationally more expensive than necessary. For example, certain types of problems, such as, e.g., monostatic scattering problems or other implicit time domain methods (where an iterative solver is used in each time step), have linearly dependent excitations—a property that conventional methods fail to exploit. Even if the excitations are linearly independent, their corresponding Krylov subspaces sometimes overlap. It is desirable to utilize this overlap to decrease the total number of iterations required to solve linear problems with multiple excitations.